% Test on Mixing Factor

%% Calling Function to calculate V(Potential), R(Charge), Err(error of consecutive steps)
N = 5001;
L = 2000;
f = 1E-4;
mf = [1E-1, 7.5E-2, 5E-2, 2.5E-2, 1E-2];
% Physical Constant Setup (In SI Unit)
a0 = 0.5291772109217;  % Bohr Radius in Angstrom 
eps = 1;               % Relative permittivity
eps_0 = 8.8541878176E-12;   % Vacuum permittivity 
A = 2.62041966677E-20;      % Area per atom
vf = 1E6; 
e = 1.6E-19; % Electronic Charge
hbar = 1.05457172647E-34;
Vmax = sqrt(pi*f*(hbar^2)*(vf^2)/(A*e^2));
Rmax = f;
Vh = 27.2113850560;


%% Function Calling
for i = 1:length(mf)
    [Va,Ra,Err] = PNGraphene(f,mf(i),N,L);
    
    % Data Conversion
    V = Va*Vh;
    R = Ra;
    
    % Data Arrangement
    N1 = length(V);
    x = linspace(-L*a0/2, L*a0/2,N1)'; 
    V = V';
    R = R';
    
    % Write to file
    filename = sprintf('pnGr(f=%1.1E,mf=%1.1E,L=%d,N=%d)',f,mf(i),L,N);
    dlmwrite([filename, 'V.txt'], [x,V]);
    dlmwrite([filename, 'R.txt'], [x,R]);
    dlmwrite([filename, 'Err.txt'], Err);
    
end

%% Data Plotting

N = 5001;
L = 2000;
f = 1E-4;
mf = [1E-1, 7.5E-2, 5E-2, 2.5E-2, 1E-2];

figure; hold on;
xlabel('$\it{Iteration (n)}$', 'interpreter', 'latex');
ylabel('$\it{|v_{n+1}-v_n|}$', 'interpreter', 'latex');

linespec={'-r', '-g', '-b', '-m','-k'};
for i = 1:length(mf)
    filename = sprintf('pnGr(f=%1.1E,mf=%1.1E,L=%d,N=%d)',f,mf(i),L,N);
    display(filename);
    Err = dlmread([filename,'Err.txt']);
    plot(Err, linespec{i});
end
legend_str = strread(num2str(mf,'mf=%1.1E\n'),'%s');
legend(legend_str);

    
    




